Complete Solution to Potential and E-Field of a sphere of radius R and a charge density ρ[r] = CC (r^2)  

Derive the electric field and electric potential both inside and outside of a sphere of radius R with a radially symmetric charge distribution given by ρ[r] = CC (r^2).  Use M as much as possible - some derivation by hand is usually needed.

					

(a)  First, determine the total charge Qinside[r] contained within the symmetrically placed sphere of radius r <= R

(*  Input cell *)
ClearAll["`*"]
ρ[r_]=CC(r^2);
Qinside[r_] = (∫_0^r)ρ[r] 4 π (rprime^2)\[DifferentialD]rprime

(b)  Find the total charge Qtotal contained in the entire sphere:

We find Qtotal by evaluating Qinside at r = R:

(*  Input cell *)
Qtotal = Qinside[R]

(c)  Determine the E field  (OverVector[EE])[r] for all r; Use the above Q functions and take advantage of the spherical symmetry.

We use Gauss’s Law.  For spherical symmetry where ρ depends only on r, Gauss’s Law simplifies to the following:

(i) Outside, the E-field at r >= R will look like that of a point charge with magnitude Qtotal.

(ii) Inside, the E-field at r <= R will be dependent only on the charge inside a sphere of radius r, i.e., Qinside[r].  This Q will look like a point charge at the origin with magnitude Qinside[r].

Both fields will be radial.   (If not convinced go to this site; it will firm up our understanding of these statements:   http://hyperphysics.phy-astr.gsu.edu/hbase/electric/elesph.html#c4).  They also discuss in detail the case where the charge is distributed uniformly throughout the sphere, i.e., ρ is a constant, i.e., proportional to (r^0) = 1).

(i) implies   OUTSIDE:   (OverVector[EE])outside[r]  = ((k Qtotal)/(r^2)) (OverHat[r])        (ii) implies   INSIDE:     (OverVector[EE])inside[r]  = ((k Qinside[r])/(r^2)) (OverHat[r])            where k = (1/(4 π (Subscript[ϵ, o])))

(d)  Enter expressions in an input cell for the magnitudes of the E fields inside and outside the sphere using Qtotal and Qinside[r];  these are functions that you found above.  Note: M remembers the Qtotal and Qinside[r] functions once the cells defining them (above) are executed.

(*  Input cell *)
EEoutside[r_] = ((k Qtotal)/(r^2))  (* magnitude only *)
EEinside[r_] = ((k Qinside[r])/(r^2))  (* magnitude only *)

(e and f) Use the above E-fields to find Voutside[r] and Vinside[r].    (I entered a bunch of assumptions to make M happy.

BTW:  I’m being sloppy about things like r > R vs. r >= R.  To M:  No Big Deal.  I usually, but not always, stick in the = sign).

(e)  First, Voutside[r] ,   r >= R

We start with V[r] - V[∞] = - (∫_∞^r)E[rprime]drprime,  where rprime is a dummy variable.  [this expression is the work per unit charge to move a test charge from infinity to the position r.  Just a reminder -- this all holds for radial dependent functions.  NOTE - sign.  If troubled by this, go here:  http://www.phys.ufl.edu/~acosta/phy2061/lectures/ElectricPotential.pdf
Therefore, to find Voutside we form the following integral (magnitude only)   Voutside[r] =  - (∫_∞^r)EEoutside[rprime] drprime = - (∫_∞^r)((k Qtotal)/(rprime^2))drprime.
rprime is a dummy variable.

(*  Input cell *)
$Assumptions =(R ≤  rprime ≤ ∞) && rprime ∈ Reals && r∈ Reals && r ≥ R;
Voutside[r_] = - (∫_∞^r)EEoutside[rprime] \[DifferentialD]rprime   (* (∫_∞^r)EEoutside[r] is defined above *)

(f)  Now find Vinside[r]   r <= R.

To find Vinside we integrate the Eoutside[r] from ∞ to R (which equals Voutside[R]) and then integrate Einside[r] from R to r; I call this contribution Vinsideonly[r]  (of course we insert the - sign).  Add the two together.  The limits on these integrals are consistent with moving the test charge from ∞ to an r inside the sphere, r < R.  Note:  I have changed the assumptions to aid the Vinsideonly[r] integration.  Recommend using Simplify near or at end.

(*  Input cell *)
$Assumptions =( 0 ≤  rprime ≤ R) && rprime ∈ Reals && ( 0 ≤  r ≤ R) && r∈ Reals;
Vinsideonly[r_] =-(∫_R^r)EEinside[rprime]\[DifferentialD]rprime
Vinside[r_] = Voutside[R] + Vinsideonly[r];
Vinside[r_] = %//Simplify

(g)  Put the two together with an If statement like this:  VV[r_] = If[ r ≥ R, Voutside[r], If[ r≤ R, Vinside[r]]], which defines the potential for all r.

(*  Input cell *)
VV[r_] = If[r ≥ R,Voutside[r],If[r≤ R,Vinside[r]]];

(h)  Do the same for E:  call it EE[r_]  (using another If statement)

EE[r_] = If[r ≥ R,EEoutside[r],If[r≤ R,EEinside[r]]];

(i)  Plot VV[r] and EE[r] together. (I used Show)

Don’t forget you have to give values for k, CC, and R for M to plot.  I used:  k = 1; CC = 1; R = 1;

(*  Input cell *)
k = 1; CC = 1; R = 1;
pp = Plot[{EE[r],VV[r]}, {r,0,6},PlotRange->{{0,6},{0,6}},PlotLegends->{"EE[r]", "V[r]"}];
line = Graphics[Line[{{R,0},{R,6}}]];
Show[pp,line,PlotRange->{{0,6},{0,6}}]  (*  The line shows the edge of the sphere, r = R *)

(j) I’d like to explore V and E for ρ[r] = CC (r^n) where we can vary n.    Reminder:  this is the n in the CHARGE DISTRIBUTION ρ[r]!

I’ve combined all the cells used above and have done the heavy lifting for you on the Manipulate which allows you to see what happens when you change n.

READ the Print Outputs carefully to understand the equations showing up in the Output.

Play with n and examine the shape of the resulting V and E functions.  A vertical line indicates where r = R.  

 If you look at the equations for Vinside, you will see that certain values of n lead to problems (try to figure out what they are).  I have tried to limit the range of r so the plots behave.
 
 Below you will be asked to do some interpretation.

(*  Execute this whole cell *)
ClearAll["`*"]
ρ[r_]=CC (r^n);
Print["ρ[r] = ",ρ[r]]
$Assumptions = 0≤ r ≤ R && r ∈ Reals && R >0 && R∈ Reals && n > -3;
Print[]
Qinside[r_] = (∫_0^r)ρ[rprime] 4 π (rprime^2)\[DifferentialD]rprime;
Print["Qinside[r] = ",Qinside[r]]

Qtotal = Qinside[R];
Print["Qtotal = Qinside[R] = ", Qtotal]
EEinside[r_] = ((k Qinside[r])/(r^2));
Print["EEinside[r] = ",EEinside[r]]
EEoutside[r_] =((k Qtotal)/(r^2)) ;
Print["EEoutside[r] = ",EEoutside[r]]
$Assumptions = r ≥ R && r ∈ Reals && R >0 && R∈ Reals && n > -3;
VVoutside[r_] = - (∫_∞^r)EEoutside[rprime]\[DifferentialD]rprime//Simplify;
Print["VVoutside[r] = ",VVoutside[r]]
VVoutside[R];
$Assumptions = 0 ≤ r ≤ R && r ∈ Reals && R >0 && R∈ Reals;
VVinside[r_] = VVoutside[R] - (∫_R^r)EEinside[rprime] \[DifferentialD]rprime;
VVinside[r_] = VVinside[r]//Simplify;
Print["VVinside[r] = ",VVinside[r]]
EE[r_] = If[r >R, EEoutside[r],If[r ≤ R,EEinside[r]]];
VVV[r_] = If[r >R, VVoutside[r],If[r ≤ R,VVinside[r]]];
k = 1; CC = 1; R = 1; rmax = 5;
Print[]
Print["Note all of these functions depend on n!"]
Print[]
Print["Here we Manipulate n to show the behavior of the magnitudes of E and V. I also show ρ[r].  The line shows r = R.  Move the slider to change n (I start with n = 5; moving the slider left to right reduces n from 5 to ~ - 1.25.  You can also hit the + and - buttons to change n by +- 0.2"]
Manipulate[ Show[Plot[{If[r>R,0,If[r≤R,(r^n)]],If[r>R,((4 π)/((3+n) (r^2))),If[r≤R,((4 π (r^(1+n)))/(3+n))]],If[r>R,((4 π)/((3+n) r)),If[r≤R,((4 π (3+n-(r^(2+n))))/((2+n) (3+n)))]]},{r,0,rmax},PlotRange-> {{0,rmax},{0,15}},PlotLegends->{"ρ[r]","EE[r]", "V[r]"}],
 Graphics[Line[{{R,0},{R,15}}]]], {n , 5,-1.25,.2,Appearance->"Open"}]

Interpret the behavior of V and E when you vary n.  Look carefully at the plots above as you vary n.

Some suggestions:  

What about V and E for r > R and all n?

Look at the behavior of E(r) with respect to the behavior of V (think E = - Grad[V]).

For r < R (inside the sphere):  Look at n = 2 (compare with above);  n = 1; n = 0 (meaning?);   0 < n < -1;  n = -1.

Anything Else?

Text Cell :  type in your interpretation(s) here.

V and E for r > R and various n:

From 	Gauss’s Law we concluded that OUTSIDE, E and therefore V looks like that of a point charge at the origin and magnitude Qtotal.  Thus, V ~ (1/r) and E ~  (1/(r^2)).  So they do not change shape with n.  (amplitude changes because of n dependence in the Q’s)

In General:  Look at the behavior of E(r) with respect to the behavior of V (think E = - Grad[V]).

For a given R, the slope of V (sign and magnitude) determines the value of E.  Note that for the n range chosen, at r = R, (dV/dr) has a maximum negative value, then continues to be negative but decreasing.  Putting in the - sign ( here E = - (dV/dr) ), E should increase with r while r < R and then decrease with r outside of the sphere (r > R).  Obviously, E has a maximum at r = R.  Although not obvious at - values of n (but very convincing at + values), the slope of V at r = 0 equals zero SO E = 0.

Look at n = 2 (compare with our first results above).   V and E look exactly the same as above (as it should).

Inside:    VVinside[r]= ((4 CC k π (-(r^(2+n))+(3+n) (R^(2+n))))/((2+n) (3+n))); n = 2 -> VVinside[r] =   ((4 CC k π (-(r^4)+ 4 (R^4)))/20).  I.e., for 0 < r < R=1, VVinside is a max at r = 0 and drops with curvature to r = R (= 1);

	EEinside[r]= ((4 CC k π (r^(1+n)))/(3+n));  n = 2 -> EEinside[r] = (4/5) CC k π (r^3)  I.e., 0 < r < R=1, EEinside is zero at r = 0 and rises as (r^3) to r = R (=1)

n = 1.  Outside r > R the curves that are consistent with V~ (1/r) and E ~ (1/(r^2)) as expected.  

Inside:    VVinside[r]= ((4 CC k π (-(r^(2+n))+(3+n) (R^(2+n))))/((2+n) (3+n))); n = 1 -> VVinside[r] =   ((4 CC k π (-(r^3)+ 4 (R^3)))/12).  I.e., for 0 < r < R=1, VVinside is a max at r = 0 and drops with curvature to r = R (= 1);

	EEinside[r]= ((4 CC k π (r^(1+n)))/(3+n));  n = 1 -> EEinside[r] = CC k π (r^2)  I.e., 0 < r < R=1, EEinside is zero at r = 0 and rises as (r^2) to r = R (=1)


 n = 0 (meaning?)  This represents a charge distribution ρ that is constant (uniform).  This leads to a linearly increasing E (perhaps you have seen this in a textbook example).  [See again:   http://hyperphysics.phy-astr.gsu.edu/hbase/electric/elesph.html#c4 ]

0 < n < -1;    Note that E is “curved” with a decreasing rate (magnitude of E with increasing r).  This can be seen by examining the r dependence of E for n < 0.  

BUT at n = -1:
Inside:    VVinside[r]= ((4 CC k π (-(r^(2+n))+(3+n) (R^(2+n))))/((2+n) (3+n))); n = -1 -> VVinside[r] =   ((4 CC k π (-r + 4 R))/12).  I.e., for 0 < r < R=1, VVinside is a max at r = 0 and drops linearly to r = R (=1)

	EEinside[r]= ((4 CC k π (r^(1+n)))/(3+n));  n = 1 -> EEinside[r] = (4/5) CC k π (r^0)  I.e., 0 < r < R=1, EEinside is CONSTANT.

-1 < n < 2;  Note the explosive behavior of E and V for certain values of n.  The most negative n I chose for the Manipulate was my attempt to keep the plot of the potential interesting.

Here is V and E for n = -.8:

Here are the functions for this particular n (-0.8):

Clear[k, CC, R]
n = -.8;  (* Reminder:  this is the n in the CHARGE DISTRIBUTION ρ[r] *)
VVinside[r]
EEinside[r]

The minus sign in the r term in VV might bother you; the R term wins.  Here is a re-plot of VVinside and EEinside (only correct for r < R) to show that all is well:

CC = 1; k = 1; R = 1; n = -.8;
EEinside[r]
VVinside[r]
Plot[{EEinside[r],VVinside[r]},{r,0,R},PlotRange-> {{0,5},{0,All}},PlotLegends->{"Einside[r]", "Vinside[r]"}]

The E and V functions only hold for r < R so the plot is terminated at r = R.  If you compare with the plot above, you see that they are the same.

For n < -2, the infinities cause so many problems that we close our eyes and walk away.  

If you have patience galore, feel free to explore. (rhyme unintended).